!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine tvd_lf (ro, mx, my, mz, bx, by, bz, en, po, gx, gy, gz, x, y, z, dx, dy, dz, dt, nx, ny, nz,                            &
   fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd,                                       &
   fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd,                                       &
   fzro_bnd, fzmx_bnd, fzmy_bnd, fzmz_bnd, fzbx_bnd, fzby_bnd, fzbz_bnd, fzen_bnd, fzpo_bnd)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz

   real(8), intent(in) :: dx, dy, dz, dt, x(nx), y(ny), z(nz)
   real(8), dimension(nx, ny, nz), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, po, gx, gy, gz
   real(8), dimension(ny, nz, n_neighs), intent(inout) ::                                                                          &
      fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd
   real(8), dimension(nx, nz, n_neighs), intent(inout) ::                                                                          &
      fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd
   real(8), dimension(ny, nz, n_neighs), intent(inout) ::                                                                          &
      fzro_bnd, fzmx_bnd, fzmy_bnd, fzmz_bnd, fzbx_bnd, fzby_bnd, fzbz_bnd, fzen_bnd, fzpo_bnd

   integer(4) :: ip1, jp1, kp1, im1, jm1, km1
   integer(4) :: sgn

   real(8) :: limiter, ch, ch2, cp2
   real(8) :: dt_dx, dt_dy, dt_dz, u1, u2, vx, vy, vz, b2, v2, c2, s2, ca2, cfx, cfy, cfz, pr, dx2, dy2, dz2

   real(8), dimension(nx, ny, nz) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre
   real(8), dimension(nx, ny, nz) :: ro_lx, mx_lx, my_lx, mz_lx, bx_lx, by_lx, bz_lx, en_lx, po_lx
   real(8), dimension(nx, ny, nz) :: ro_rx, mx_rx, my_rx, mz_rx, bx_rx, by_rx, bz_rx, en_rx, po_rx
   real(8), dimension(nx, ny, nz) :: ro_ly, mx_ly, my_ly, mz_ly, bx_ly, by_ly, bz_ly, en_ly, po_ly
   real(8), dimension(nx, ny, nz) :: ro_ry, mx_ry, my_ry, mz_ry, bx_ry, by_ry, bz_ry, en_ry, po_ry
   real(8), dimension(nx, ny, nz) :: ro_lz, mx_lz, my_lz, mz_lz, bx_lz, by_lz, bz_lz, en_lz, po_lz
   real(8), dimension(nx, ny, nz) :: ro_rz, mx_rz, my_rz, mz_rz, bx_rz, by_rz, bz_rz, en_rz, po_rz
   real(8), dimension(nx, ny, nz) :: fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo
   real(8), dimension(nx, ny, nz) :: fyro, fymx, fymy, fymz, fybx, fyby, fybz, fyen, fypo
   real(8), dimension(nx, ny, nz) :: fzro, fzmx, fzmy, fzmz, fzbx, fzby, fzbz, fzen, fzpo
   real(8), dimension(nx, ny, nz) :: bx_resis, by_resis, bz_resis, en_resis

   real(8) :: fxro_l, fxmx_l, fxmy_l, fxmz_l, fxbx_l, fxby_l, fxbz_l, fxen_l, fxpo_l
   real(8) :: fxro_r, fxmx_r, fxmy_r, fxmz_r, fxbx_r, fxby_r, fxbz_r, fxen_r, fxpo_r
   real(8) :: fyro_l, fymx_l, fymy_l, fymz_l, fybx_l, fyby_l, fybz_l, fyen_l, fypo_l
   real(8) :: fyro_r, fymx_r, fymy_r, fymz_r, fybx_r, fyby_r, fybz_r, fyen_r, fypo_r
   real(8) :: fzro_l, fzmx_l, fzmy_l, fzmz_l, fzbx_l, fzby_l, fzbz_l, fzen_l, fzpo_l
   real(8) :: fzro_r, fzmx_r, fzmy_r, fzmz_r, fzbx_r, fzby_r, fzbz_r, fzen_r, fzpo_r

!----------------------------------------------------------------------------------------------------------------------------------|
!  left and right reconstructed value
!----------------------------------------------------------------------------------------------------------------------------------|
   dt_dx = dt / dx
   dt_dy = dt / dy
   dt_dz = dt / dz
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0
   dz2 = dz * 2.0d0
   ch = cfl * min(dx, dy, dz) / dt * div_clean_flag
   ch2 = ch * ch
!   cp2 = ch * cd / 6.0d0 + eps
   cp2 = -dt * ch2 / log (cd) + eps

   do k = 2, nz - 1
      kp1 = k + 1
      km1 = k - 1
      do j = 2, ny - 1
         jp1 = j + 1
         jm1 = j - 1
         do i = 2, nx - 1
            ip1 = i + 1
            im1 = i - 1

            u1 = ro(i, j, k) - ro(im1, j, k)
            u2 = ro(ip1, j, k) - ro(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            ro_lx(i, j, k) = ro(i, j, k) + 0.5d0 * limiter
            ro_rx(i, j, k) = ro(i, j, k) - 0.5d0 * limiter

            u1 = mx(i, j, k) - mx(im1, j, k)
            u2 = mx(ip1, j, k) - mx(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            mx_lx(i, j, k) = mx(i, j, k) + 0.5d0 * limiter
            mx_rx(i, j, k) = mx(i, j, k) - 0.5d0 * limiter

            u1 = my(i, j, k) - my(im1, j, k)
            u2 = my(ip1, j, k) - my(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            my_lx(i, j, k) = my(i, j, k) + 0.5d0 * limiter
            my_rx(i, j, k) = my(i, j, k) - 0.5d0 * limiter

            u1 = mz(i, j, k) - mz(im1, j, k)
            u2 = mz(ip1, j, k) - mz(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            mz_lx(i, j, k) = mz(i, j, k) + 0.5d0 * limiter
            mz_rx(i, j, k) = mz(i, j, k) - 0.5d0 * limiter

            u1 = bx(i, j, k) - bx(im1, j, k)
            u2 = bx(ip1, j, k) - bx(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            bx_lx(i, j, k) = bx(i, j, k) + 0.5d0 * limiter
            bx_rx(i, j, k) = bx(i, j, k) - 0.5d0 * limiter

            u1 = by(i, j, k) - by(im1, j, k)
            u2 = by(ip1, j, k) - by(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            by_lx(i, j, k) = by(i, j, k) + 0.5d0 * limiter
            by_rx(i, j, k) = by(i, j, k) - 0.5d0 * limiter

            u1 = bz(i, j, k) - bz(im1, j, k)
            u2 = bz(ip1, j, k) - bz(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            bz_lx(i, j, k) = bz(i, j, k) + 0.5d0 * limiter
            bz_rx(i, j, k) = bz(i, j, k) - 0.5d0 * limiter

            u1 = en(i, j, k) - en(im1, j, k)
            u2 = en(ip1, j, k) - en(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            en_lx(i, j, k) = en(i, j, k) + 0.5d0 * limiter
            en_rx(i, j, k) = en(i, j, k) - 0.5d0 * limiter

            u1 = po(i, j, k) - po(im1, j, k)
            u2 = po(ip1, j, k) - po(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            po_lx(i, j, k) = po(i, j, k) + 0.5d0 * limiter
            po_rx(i, j, k) = po(i, j, k) - 0.5d0 * limiter

            u1 = ro(i, j, k) - ro(i, jm1, k)
            u2 = ro(i, jp1, k) - ro(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            ro_ly(i, j, k) = ro(i, j, k) + 0.5d0 * limiter
            ro_ry(i, j, k) = ro(i, j, k) - 0.5d0 * limiter

            u1 = mx(i, j, k) - mx(i, jm1, k)
            u2 = mx(i, jp1, k) - mx(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            mx_ly(i, j, k) = mx(i, j, k) + 0.5d0 * limiter
            mx_ry(i, j, k) = mx(i, j, k) - 0.5d0 * limiter

            u1 = my(i, j, k) - my(i, jm1, k)
            u2 = my(i, jp1, k) - my(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            my_ly(i, j, k) = my(i, j, k) + 0.5d0 * limiter
            my_ry(i, j, k) = my(i, j, k) - 0.5d0 * limiter

            u1 = mz(i, j, k) - mz(i, jm1, k)
            u2 = mz(i, jp1, k) - mz(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            mz_ly(i, j, k) = mz(i, j, k) + 0.5d0 * limiter
            mz_ry(i, j, k) = mz(i, j, k) - 0.5d0 * limiter

            u1 = bx(i, j, k) - bx(i, jm1, k)
            u2 = bx(i, jp1, k) - bx(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            bx_ly(i, j, k) = bx(i, j, k) + 0.5d0 * limiter
            bx_ry(i, j, k) = bx(i, j, k) - 0.5d0 * limiter

            u1 = by(i, j, k) - by(i, jm1, k)
            u2 = by(i, jp1, k) - by(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            by_ly(i, j, k) = by(i, j, k) + 0.5d0 * limiter
            by_ry(i, j, k) = by(i, j, k) - 0.5d0 * limiter

            u1 = bz(i, j, k) - bz(i, jm1, k)
            u2 = bz(i, jp1, k) - bz(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            bz_ly(i, j, k) = bz(i, j, k) + 0.5d0 * limiter
            bz_ry(i, j, k) = bz(i, j, k) - 0.5d0 * limiter

            u1 = en(i, j, k) - en(i, jm1, k)
            u2 = en(i, jp1, k) - en(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            en_ly(i, j, k) = en(i, j, k) + 0.5d0 * limiter
            en_ry(i, j, k) = en(i, j, k) - 0.5d0 * limiter

            u1 = po(i, j, k) - po(i, jm1, k)
            u2 = po(i, jp1, k) - po(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            po_ly(i, j, k) = po(i, j, k) + 0.5d0 * limiter
            po_ry(i, j, k) = po(i, j, k) - 0.5d0 * limiter

            u1 = ro(i, j, k) - ro(i, j, km1)
            u2 = ro(i, j, kp1) - ro(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            ro_lz(i, j, k) = ro(i, j, k) + 0.5d0 * limiter
            ro_rz(i, j, k) = ro(i, j, k) - 0.5d0 * limiter

            u1 = mx(i, j, k) - mx(i, j, km1)
            u2 = mx(i, j, kp1) - mx(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            mx_lz(i, j, k) = mx(i, j, k) + 0.5d0 * limiter
            mx_rz(i, j, k) = mx(i, j, k) - 0.5d0 * limiter

            u1 = my(i, j, k) - my(i, j, km1)
            u2 = my(i, j, kp1) - my(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            my_lz(i, j, k) = my(i, j, k) + 0.5d0 * limiter
            my_rz(i, j, k) = my(i, j, k) - 0.5d0 * limiter

            u1 = mz(i, j, k) - mz(i, j, km1)
            u2 = mz(i, j, kp1) - mz(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            mz_lz(i, j, k) = mz(i, j, k) + 0.5d0 * limiter
            mz_rz(i, j, k) = mz(i, j, k) - 0.5d0 * limiter

            u1 = bx(i, j, k) - bx(i, j, km1)
            u2 = bx(i, j, kp1) - bx(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            bx_lz(i, j, k) = bx(i, j, k) + 0.5d0 * limiter
            bx_rz(i, j, k) = bx(i, j, k) - 0.5d0 * limiter

            u1 = by(i, j, k) - by(i, j, km1)
            u2 = by(i, j, kp1) - by(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            by_lz(i, j, k) = by(i, j, k) + 0.5d0 * limiter
            by_rz(i, j, k) = by(i, j, k) - 0.5d0 * limiter

            u1 = bz(i, j, k) - bz(i, j, km1)
            u2 = bz(i, j, kp1) - bz(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            bz_lz(i, j, k) = bz(i, j, k) + 0.5d0 * limiter
            bz_rz(i, j, k) = bz(i, j, k) - 0.5d0 * limiter

            u1 = en(i, j, k) - en(i, j, km1)
            u2 = en(i, j, kp1) - en(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            en_lz(i, j, k) = en(i, j, k) + 0.5d0 * limiter
            en_rz(i, j, k) = en(i, j, k) - 0.5d0 * limiter

            u1 = po(i, j, k) - po(i, j, km1)
            u2 = po(i, j, kp1) - po(i, j, k)
            if (u1 .ge. 0.0d0) then
               sgn = 1
            else
               sgn = -1
            endif
            limiter = sgn * max(0.0d0, min(abs(u1), sgn * u2))
            po_lz(i, j, k) = po(i, j, k) + 0.5d0 * limiter
            po_rz(i, j, k) = po(i, j, k) - 0.5d0 * limiter
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  predictor-corrector step
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx
      vx = mx_lx(i, j, k) / ro_lx(i, j, k)
      vy = my_lx(i, j, k) / ro_lx(i, j, k)
      vz = mz_lx(i, j, k) / ro_lx(i, j, k)
      b2 = bx_lx(i, j, k) * bx_lx(i, j, k) + by_lx(i, j, k) * by_lx(i, j, k) + bz_lx(i, j, k) * bz_lx(i, j, k)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_lx(i, j, k) - v2 * ro_lx(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
      fxro_l = mx_lx(i, j, k)
      fxmx_l = vx * mx_lx(i, j, k) + pr + b2 / 2.0d0 - bx_lx(i, j, k) * bx_lx(i, j, k)
      fxmy_l = vx * my_lx(i, j, k) - bx_lx(i, j, k) * by_lx(i, j, k)
      fxmz_l = vx * mz_lx(i, j, k) - bx_lx(i, j, k) * bz_lx(i, j, k)
      fxbx_l = po_lx(i, j, k)
      fxby_l = vx * by_lx(i, j, k) - bx_lx(i, j, k) * vy
      fxbz_l = vx * bz_lx(i, j, k) - bx_lx(i, j, k) * vz
      fxen_l = (en_lx(i, j, k) + pr + b2 / 2.0d0) * vx -                                                                           &
         (bx_lx(i, j, k) * vx + by_lx(i, j, k) * vy + bz_lx(i, j, k) * vz) * bx_lx(i, j, k)
      fxpo_l = ch2 * bx_lx(i, j, k)

      vx = mx_rx(i, j, k) / ro_rx(i, j, k)
      vy = my_rx(i, j, k) / ro_rx(i, j, k)
      vz = mz_rx(i, j, k) / ro_rx(i, j, k)
      b2 = bx_rx(i, j, k) * bx_rx(i, j, k) + by_rx(i, j, k) * by_rx(i, j, k) + bz_rx(i, j, k) * bz_rx(i, j, k)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_rx(i, j, k) - v2 * ro_rx(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
      fxro_r = mx_rx(i, j, k)
      fxmx_r = vx * mx_rx(i, j, k) + pr + b2 / 2.0d0 - bx_rx(i, j, k) * bx_rx(i, j, k)
      fxmy_r = vx * my_rx(i, j, k) - bx_rx(i, j, k) * by_rx(i, j, k)
      fxmz_r = vx * mz_rx(i, j, k) - bx_rx(i, j, k) * bz_rx(i, j, k)
      fxbx_r = po_rx(i, j, k)
      fxby_r = vx * by_rx(i, j, k) - bx_rx(i, j, k) * vy
      fxbz_r = vx * bz_rx(i, j, k) - bx_rx(i, j, k) * vz
      fxen_r = (en_rx(i, j, k) + pr + b2 / 2.0d0) * vx -                                                                           &
         (bx_rx(i, j, k) * vx + by_rx(i, j, k) * vy + bz_rx(i, j, k) * vz) * bx_rx(i, j, k)
      fxpo_r = ch2 * bx_rx(i, j, k)

      vx = mx_ly(i, j, k) / ro_ly(i, j, k)
      vy = my_ly(i, j, k) / ro_ly(i, j, k)
      vz = mz_ly(i, j, k) / ro_ly(i, j, k)
      b2 = bx_ly(i, j, k) * bx_ly(i, j, k) + by_ly(i, j, k) * by_ly(i, j, k) + bz_ly(i, j, k) * bz_ly(i, j, k)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_ly(i, j, k) - v2 * ro_ly(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
      fyro_l = my_ly(i, j, k)
      fymx_l = vy * mx_ly(i, j, k) - by_ly(i, j, k) * bx_ly(i, j, k)
      fymy_l = vy * my_ly(i, j, k) + pr + b2 / 2.0d0 - by_ly(i, j, k) * by_ly(i, j, k)
      fymz_l = vy * mz_ly(i, j, k) - by_ly(i, j, k) * bz_ly(i, j, k)
      fybx_l = vy * bx_ly(i, j, k) - by_ly(i, j, k) * vx
      fyby_l = po_ly(i, j, k)
      fybz_l = vy * bz_ly(i, j, k) - by_ly(i, j, k) * vz
      fyen_l = (en_ly(i, j, k) + pr + b2 / 2.0d0) * vy -                                                                           &
         (bx_ly(i, j, k) * vx + by_ly(i, j, k) * vy + bz_ly(i, j, k) * vz) * by_ly(i, j, k)
      fypo_l = ch2 * by_ly(i, j, k)

      vx = mx_ry(i, j, k) / ro_ry(i, j, k)
      vy = my_ry(i, j, k) / ro_ry(i, j, k)
      vz = mz_ry(i, j, k) / ro_ry(i, j, k)
      b2 = bx_ry(i, j, k) * bx_ry(i, j, k) + by_ry(i, j, k) * by_ry(i, j, k) + bz_ry(i, j, k) * bz_ry(i, j, k)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_ry(i, j, k) - v2 * ro_ry(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
      fyro_r = my_ry(i, j, k)
      fymx_r = vy * mx_ry(i, j, k) - by_ry(i, j, k) * bx_ry(i, j, k)
      fymy_r = vy * my_ry(i, j, k) + pr + b2 / 2.0d0 - by_ry(i, j, k) * by_ry(i, j, k)
      fymz_r = vy * mz_ry(i, j, k) - by_ry(i, j, k) * bz_ry(i, j, k)
      fybx_r = vy * bx_ry(i, j, k) - by_ry(i, j, k) * vx
      fyby_r = po_ry(i, j, k)
      fybz_r = vy * bz_ry(i, j, k) - by_ry(i, j, k) * vz
      fyen_r = (en_ry(i, j, k) + pr + b2 / 2.0d0) * vy -                                                                           &
         (bx_ry(i, j, k) * vx + by_ry(i, j, k) * vy + bz_ry(i, j, k) * vz) * by_ry(i, j, k)
      fypo_r = ch2 * by_ry(i, j, k)

      vx = mx_lz(i, j, k) / ro_lz(i, j, k)
      vy = my_lz(i, j, k) / ro_lz(i, j, k)
      vz = mz_lz(i, j, k) / ro_lz(i, j, k)
      b2 = bx_lz(i, j, k) * bx_lz(i, j, k) + by_lz(i, j, k) * by_lz(i, j, k) + bz_lz(i, j, k) * bz_lz(i, j, k)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_lz(i, j, k) - v2 * ro_lz(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
      fzro_l = mz_lz(i, j, k)
      fzmx_l = vz * mx_lz(i, j, k) - bz_lz(i, j, k) * bx_lz(i, j, k)
      fzmy_l = vz * my_lz(i, j, k) - bz_lz(i, j, k) * by_lz(i, j, k)
      fzmz_l = vz * mz_lz(i, j, k) + pr + b2 / 2.0d0 - bz_lz(i, j, k) * bz_lz(i, j, k)
      fzbx_l = vz * bx_lz(i, j, k) - bz_lz(i, j, k) * vx
      fzby_l = vz * by_lz(i, j, k) - bz_lz(i, j, k) * vy
      fzbz_l = po_lz(i, j, k)
      fzen_l = (en_lz(i, j, k) + pr + b2 / 2.0d0) * vz -                                                                           &
         (bx_lz(i, j, k) * vx + by_lz(i, j, k) * vy + bz_lz(i, j, k) * vz) * bz_lz(i, j, k)
      fzpo_l = ch2 * bz_lz(i, j, k)

      vx = mx_rz(i, j, k) / ro_rz(i, j, k)
      vy = my_rz(i, j, k) / ro_rz(i, j, k)
      vz = mz_rz(i, j, k) / ro_rz(i, j, k)
      b2 = bx_rz(i, j, k) * bx_rz(i, j, k) + by_rz(i, j, k) * by_rz(i, j, k) + bz_rz(i, j, k) * bz_rz(i, j, k)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_rz(i, j, k) - v2 * ro_rz(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
      fzro_l = mz_rz(i, j, k)
      fzmx_l = vz * mx_rz(i, j, k) - bz_rz(i, j, k) * bx_rz(i, j, k)
      fzmy_l = vz * my_rz(i, j, k) - bz_rz(i, j, k) * by_rz(i, j, k)
      fzmz_l = vz * mz_rz(i, j, k) + pr + b2 / 2.0d0 - bz_rz(i, j, k) * bz_rz(i, j, k)
      fzbx_l = vz * bx_rz(i, j, k) - bz_rz(i, j, k) * vx
      fzby_l = vz * by_rz(i, j, k) - bz_rz(i, j, k) * vy
      fzbz_l = po_rz(i, j, k)
      fzen_l = (en_rz(i, j, k) + pr + b2 / 2.0d0) * vz -                                                                           &
         (bx_rz(i, j, k) * vx + by_rz(i, j, k) * vy + bz_rz(i, j, k) * vz) * bz_rz(i, j, k)
      fzpo_l = ch2 * bz_rz(i, j, k)

      ro_pre(i, j, k) = ro(i, j, k) - 0.5d0 * dt_dx * (fxro_l - fxro_r) -                                                          &
                                      0.5d0 * dt_dy * (fyro_l - fyro_r) -                                                          &
                                      0.5d0 * dt_dz * (fzro_l - fzro_r)
      mx_pre(i, j, k) = mx(i, j, k) - 0.5d0 * dt_dx * (fxmx_l - fxmx_r) -                                                          &
                                      0.5d0 * dt_dy * (fymx_l - fymx_r) -                                                          &
                                      0.5d0 * dt_dz * (fzmx_l - fzmx_r)
      my_pre(i, j, k) = my(i, j, k) - 0.5d0 * dt_dx * (fxmy_l - fxmy_r) -                                                          &
                                      0.5d0 * dt_dy * (fymy_l - fymy_r) -                                                          &
                                      0.5d0 * dt_dz * (fzmy_l - fzmy_r)
      mz_pre(i, j, k) = mz(i, j, k) - 0.5d0 * dt_dx * (fxmz_l - fxmz_r) -                                                          &
                                      0.5d0 * dt_dy * (fymz_l - fymz_r) -                                                          &
                                      0.5d0 * dt_dz * (fzmz_l - fzmz_r)
      bx_pre(i, j, k) = bx(i, j, k) - 0.5d0 * dt_dx * (fxbx_l - fxbx_r) -                                                          &
                                      0.5d0 * dt_dy * (fybx_l - fybx_r) -                                                          &
                                      0.5d0 * dt_dz * (fzbx_l - fzbx_r)
      by_pre(i, j, k) = by(i, j, k) - 0.5d0 * dt_dx * (fxby_l - fxby_r) -                                                          &
                                      0.5d0 * dt_dy * (fyby_l - fyby_r) -                                                          &
                                      0.5d0 * dt_dz * (fzby_l - fzby_r)
      bz_pre(i, j, k) = bz(i, j, k) - 0.5d0 * dt_dx * (fxbz_l - fxbz_r) -                                                          &
                                      0.5d0 * dt_dy * (fybz_l - fybz_r) -                                                          &
                                      0.5d0 * dt_dz * (fzbz_l - fzbz_r)
      en_pre(i, j, k) = en(i, j, k) - 0.5d0 * dt_dx * (fxen_l - fxen_r) -                                                          &
                                      0.5d0 * dt_dy * (fyen_l - fyen_r) -                                                          &
                                      0.5d0 * dt_dz * (fzen_l - fzen_r)
      po_pre(i, j, k) = po(i, j, k) - 0.5d0 * dt_dx * (fxpo_l - fxpo_r) -                                                          &
                                      0.5d0 * dt_dy * (fypo_l - fypo_r) -                                                          &
                                      0.5d0 * dt_dz * (fzpo_l - fzpo_r)
   enddo
   enddo
   enddo

   call resistivity(bx_resis, by_resis, bz_resis, en_resis, ro, mx, my, mz, bx, by, bz, en, x, y, z, dx, dy, dz, nx, ny, nz)

   call source_term (mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre,                                               &
      ro, mx, my, mz, bx, by, bz, po, gx, gy, gz, bx_resis, by_resis, bz_resis, en_resis,                                          &
      ch2, cp2, dx, dy, dz, dt / 2.0d0, nx, ny, nz)

!----------------------------------------------------------------------------------------------------------------------------------|
!  left and righe reconstructed value
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx

      limiter = ro_pre(i, j, k) - ro(i, j, k)
      ro_lx(i, j, k) = ro_lx(i, j, k) + limiter
      ro_rx(i, j, k) = ro_rx(i, j, k) + limiter
      ro_ly(i, j, k) = ro_ly(i, j, k) + limiter
      ro_ry(i, j, k) = ro_ry(i, j, k) + limiter
      ro_lz(i, j, k) = ro_lz(i, j, k) + limiter
      ro_rz(i, j, k) = ro_rz(i, j, k) + limiter

      limiter = mx_pre(i, j, k) - mx(i, j, k)
      mx_lx(i, j, k) = mx_lx(i, j, k) + limiter
      mx_rx(i, j, k) = mx_rx(i, j, k) + limiter
      mx_ly(i, j, k) = mx_ly(i, j, k) + limiter
      mx_ry(i, j, k) = mx_ry(i, j, k) + limiter
      mx_lz(i, j, k) = mx_lz(i, j, k) + limiter
      mx_rz(i, j, k) = mx_rz(i, j, k) + limiter

      limiter = my_pre(i, j, k) - my(i, j, k)
      my_lx(i, j, k) = my_lx(i, j, k) + limiter
      my_rx(i, j, k) = my_rx(i, j, k) + limiter
      my_ly(i, j, k) = my_ly(i, j, k) + limiter
      my_ry(i, j, k) = my_ry(i, j, k) + limiter
      my_lz(i, j, k) = my_lz(i, j, k) + limiter
      my_rz(i, j, k) = my_rz(i, j, k) + limiter

      limiter = mz_pre(i, j, k) - mz(i, j, k)
      mz_lx(i, j, k) = mz_lx(i, j, k) + limiter
      mz_rx(i, j, k) = mz_rx(i, j, k) + limiter
      mz_ly(i, j, k) = mz_ly(i, j, k) + limiter
      mz_ry(i, j, k) = mz_ry(i, j, k) + limiter
      mz_lz(i, j, k) = mz_lz(i, j, k) + limiter
      mz_rz(i, j, k) = mz_rz(i, j, k) + limiter

      limiter = bx_pre(i, j, k) - bx(i, j, k)
      bx_lx(i, j, k) = bx_lx(i, j, k) + limiter
      bx_rx(i, j, k) = bx_rx(i, j, k) + limiter
      bx_ly(i, j, k) = bx_ly(i, j, k) + limiter
      bx_ry(i, j, k) = bx_ry(i, j, k) + limiter
      bx_lz(i, j, k) = bx_lz(i, j, k) + limiter
      bx_rz(i, j, k) = bx_rz(i, j, k) + limiter

      limiter = by_pre(i, j, k) - by(i, j, k)
      by_lx(i, j, k) = by_lx(i, j, k) + limiter
      by_rx(i, j, k) = by_rx(i, j, k) + limiter
      by_ly(i, j, k) = by_ly(i, j, k) + limiter
      by_ry(i, j, k) = by_ry(i, j, k) + limiter
      by_lz(i, j, k) = by_lz(i, j, k) + limiter
      by_rz(i, j, k) = by_rz(i, j, k) + limiter

      limiter = bz_pre(i, j, k) - bz(i, j, k)
      bz_lx(i, j, k) = bz_lx(i, j, k) + limiter
      bz_rx(i, j, k) = bz_rx(i, j, k) + limiter
      bz_ly(i, j, k) = bz_ly(i, j, k) + limiter
      bz_ry(i, j, k) = bz_ry(i, j, k) + limiter
      bz_lz(i, j, k) = bz_lz(i, j, k) + limiter
      bz_rz(i, j, k) = bz_rz(i, j, k) + limiter

      limiter = en_pre(i, j, k) - en(i, j, k)
      en_lx(i, j, k) = en_lx(i, j, k) + limiter
      en_rx(i, j, k) = en_rx(i, j, k) + limiter
      en_ly(i, j, k) = en_ly(i, j, k) + limiter
      en_ry(i, j, k) = en_ry(i, j, k) + limiter
      en_lz(i, j, k) = en_lz(i, j, k) + limiter
      en_lz(i, j, k) = en_rz(i, j, k) + limiter

      limiter = po_pre(i, j, k) - po(i, j, k)
      po_lx(i, j, k) = po_lx(i, j, k) + limiter
      po_rx(i, j, k) = po_rx(i, j, k) + limiter
      po_ly(i, j, k) = po_ly(i, j, k) + limiter
      po_ry(i, j, k) = po_ry(i, j, k) + limiter
      po_lz(i, j, k) = po_lz(i, j, k) + limiter
      po_rz(i, j, k) = po_rz(i, j, k) + limiter

   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  Riemann solver
!----------------------------------------------------------------------------------------------------------------------------------|
   if (riemann_solver_flag .eq. 1) then
      call hlld_x (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen,                                                                 &
         ro_lx, mx_lx, my_lx, mz_lx, bx_lx, by_lx, bz_lx, en_lx, po_lx,                                                            &
         ro_rx, mx_rx, my_rx, mz_rx, bx_rx, by_rx, bz_rx, en_rx, po_rx, nx, ny, nz)
      call hlld_y (fyro, fymy, fymx, fymz, fyby, fybx, fybz, fyen,                                                                 &
         ro_ly, my_ly, mx_ly, mz_ly, by_ly, bx_ly, bz_ly, en_ly, po_ly,                                                            &
         ro_ry, my_ry, mx_ry, mz_ry, by_ry, bx_ry, bz_ry, en_ry, po_ry, nx, ny, nz)
      call hlld_z (fzro, fzmz, fzmy, fzmx, fzbz, fzby, fzbx, fzen,                                                                 &
         ro_lz, mz_lz, my_lz, mx_lz, bz_lz, by_lz, bx_lz, en_lz, po_lz,                                                            &
         ro_rz, mz_rz, my_rz, mx_rz, bz_rz, by_rz, bx_rz, en_rz, po_rz, nx, ny, nz)
   elseif (riemann_solver_flag .eq. 2) then
      call hllc_x (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen,                                                                 &
         ro_lx, mx_lx, my_lx, mz_lx, bx_lx, by_lx, bz_lx, en_lx, po_lx,                                                            &
         ro_rx, mx_rx, my_rx, mz_rx, bx_rx, by_rx, bz_rx, en_rx, po_rx, nx, ny, nz)
      call hllc_y (fyro, fymy, fymx, fymz, fyby, fybx, fybz, fyen,                                                                 &
         ro_ly, my_ly, mx_ly, mz_ly, by_ly, bx_ly, bz_ly, en_ly, po_ly,                                                            &
         ro_ry, my_ry, mx_ry, mz_ry, by_ry, bx_ry, bz_ry, en_ry, po_ry, nx, ny, nz)
      call hllc_z (fzro, fzmz, fzmy, fzmx, fzbz, fzby, fzbx, fzen,                                                                 &
         ro_lz, mz_lz, my_lz, mx_lz, bz_lz, by_lz, bx_lz, en_lz, po_lz,                                                            &
         ro_rz, mz_rz, my_rz, mx_rz, bz_rz, by_rz, bx_rz, en_rz, po_rz, nx, ny, nz)
   elseif (riemann_solver_flag .eq. 3) then
      call roe_x  (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen,                                                                 &
         ro_lx, mx_lx, my_lx, mz_lx, bx_lx, by_lx, bz_lx, en_lx, po_lx,                                                            &
         ro_rx, mx_rx, my_rx, mz_rx, bx_rx, by_rx, bz_rx, en_rx, po_rx, nx, ny, nz)
      call roe_y  (fyro, fymy, fymx, fymz, fyby, fybx, fybz, fyen,                                                                 &
         ro_ly, my_ly, mx_ly, mz_ly, by_ly, bx_ly, bz_ly, en_ly, po_ly,                                                            &
         ro_ry, my_ry, mx_ry, mz_ry, by_ry, bx_ry, bz_ry, en_ry, po_ry, nx, ny, nz)
      call roe_z  (fzro, fzmz, fzmy, fzmx, fzbz, fzby, fzbx, fzen,                                                                 &
         ro_lz, mz_lz, my_lz, mx_lz, bz_lz, by_lz, bx_lz, en_lz, po_lz,                                                            &
         ro_rz, mz_rz, my_rz, mx_rz, bz_rz, by_rz, bx_rz, en_rz, po_rz, nx, ny, nz)
   endif

   do k = 1, nz - 1
      kp1 = k + 1
      do j = 1, ny - 1
         jp1 = j + 1
         do i = 1, nx - 1
            ip1 = i + 1

            vx = mx_lx(i, j, k) / ro_lx(i, j, k)
            vy = my_lx(i, j, k) / ro_lx(i, j, k)
            vz = mz_lx(i, j, k) / ro_lx(i, j, k)
            b2 = bx_lx(i, j, k) * bx_lx(i, j, k) + by_lx(i, j, k) * by_lx(i, j, k) + bz_lx(i, j, k) * bz_lx(i, j, k)
            v2 = vx * vx + vy * vy + vz * vz
            pr = (en_lx(i, j, k) - v2 * ro_lx(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
            c2 = gm * pr
            s2 = c2 + b2
            ca2 = bx_lx(i, j, k) * bx_lx(i, j, k)
            cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_lx(i, j, k) / 2.0d0)
            max_speed = cfx + abs(vx)

            if (riemann_solver_flag .eq. 0) then
               fxro_l = mx_lx(i, j, k)
               fxmx_l = vx * mx_lx(i, j, k) + pr + b2 / 2.0d0 - bx_lx(i, j, k) * bx_lx(i, j, k)
               fxmy_l = vx * my_lx(i, j, k) - bx_lx(i, j, k) * by_lx(i, j, k)
               fxmz_l = vx * mz_lx(i, j, k) - bx_lx(i, j, k) * bz_lx(i, j, k)
               fxbx_l = po_lx(i, j, k)
               fxby_l = vx * by_lx(i, j, k) - bx_lx(i, j, k) * vy
               fxbz_l = vx * bz_lx(i, j, k) - bx_lx(i, j, k) * vz
               fxen_l = (en_lx(i, j, k) + pr + b2 / 2.0d0) * vx -                                                                  &
                  (bx_lx(i, j, k) * vx + by_lx(i, j, k) * vy + bz_lx(i, j, k) * vz) * bx_lx(i, j, k)
            endif
            fxpo_l = ch2 * bx_lx(i, j, k)

            vx = mx_rx(ip1, j, k) / ro_rx(ip1, j, k)
            vy = my_rx(ip1, j, k) / ro_rx(ip1, j, k)
            vz = mz_rx(ip1, j, k) / ro_rx(ip1, j, k)
            b2 = bx_rx(ip1, j, k) * bx_rx(ip1, j, k) + by_rx(ip1, j, k) * by_rx(ip1, j, k) + bz_rx(ip1, j, k) * bz_rx(ip1, j, k)
            v2 = vx * vx + vy * vy + vz * vz
            pr = (en_rx(ip1, j, k) - v2 * ro_rx(ip1, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
            c2 = gm * pr
            s2 = c2 + b2
            ca2 = bx_rx(ip1, j, k) * bx_rx(ip1, j, k)
            cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_rx(ip1, j, k) / 2.0d0)
            max_speed = max(cfx + abs(vx), max_speed)

            if (riemann_solver_flag .eq. 0) then
               fxro_r = mx_rx(ip1, j, k)
               fxmx_r = vx * mx_rx(ip1, j, k) + pr + b2 / 2.0d0 - bx_rx(ip1, j, k) * bx_rx(ip1, j, k)
               fxmy_r = vx * my_rx(ip1, j, k) - bx_rx(ip1, j, k) * by_rx(ip1, j, k)
               fxmz_r = vx * mz_rx(ip1, j, k) - bx_rx(ip1, j, k) * bz_rx(ip1, j, k)
               fxbx_r = po_rx(i, j, k)
               fxby_r = vx * by_rx(ip1, j, k) - bx_rx(ip1, j, k) * vy
               fxbz_r = vx * bz_rx(ip1, j, k) - bx_rx(ip1, j, k) * vz
               fxen_r = (en_rx(ip1, j, k) + pr + b2 / 2.0d0) * vx -                                                                &
                  (bx_rx(ip1, j, k) * vx + by_rx(ip1, j, k) * vy + bz_rx(ip1, j, k) * vz) * bx_rx(ip1, j, k)
            endif
            fxpo_r = ch2 * bx_rx(ip1, j, k)
            
            if (riemann_solver_flag .eq. 0) then
               fxro(i, j, k) = 0.5d0 * (fxro_l + fxro_r - max_speed * (ro_rx(ip1, j, k) - ro_lx(i, j, k)))
               fxmx(i, j, k) = 0.5d0 * (fxmx_l + fxmx_r - max_speed * (mx_rx(ip1, j, k) - mx_lx(i, j, k)))
               fxmy(i, j, k) = 0.5d0 * (fxmy_l + fxmy_r - max_speed * (my_rx(ip1, j, k) - my_lx(i, j, k)))
               fxmz(i, j, k) = 0.5d0 * (fxmz_l + fxmz_r - max_speed * (mz_rx(ip1, j, k) - mz_lx(i, j, k)))
               fxbx(i, j, k) = 0.5d0 * (fxbx_l + fxbx_r - max_speed * (bx_rx(ip1, j, k) - bx_lx(i, j, k)))
               fxby(i, j, k) = 0.5d0 * (fxby_l + fxby_r - max_speed * (by_rx(ip1, j, k) - by_lx(i, j, k)))
               fxbz(i, j, k) = 0.5d0 * (fxbz_l + fxbz_r - max_speed * (bz_rx(ip1, j, k) - bz_lx(i, j, k)))
               fxen(i, j, k) = 0.5d0 * (fxen_l + fxen_r - max_speed * (en_rx(ip1, j, k) - en_lx(i, j, k)))
            endif
            fxpo(i, j, k) = 0.5d0 * (fxpo_l + fxpo_r - max_speed * (po_rx(ip1, j, k) - po_lx(i, j, k)))

            vx = mx_ly(i, j, k) / ro_ly(i, j, k)
            vy = my_ly(i, j, k) / ro_ly(i, j, k)
            vz = mz_ly(i, j, k) / ro_ly(i, j, k)
            b2 = bx_ly(i, j, k) * bx_ly(i, j, k) + by_ly(i, j, k) * by_ly(i, j, k) + bz_ly(i, j, k) * bz_ly(i, j, k)
            v2 = vx * vx + vy * vy + vz * vz
            pr = (en_ly(i, j, k) - v2 * ro_ly(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
            c2 = gm * pr
            s2 = c2 + b2
            ca2 = by_ly(i, j, k) * by_ly(i, j, k)
            cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_ly(i, j, k) / 2.0d0)
            max_speed = cfy + abs(vy)

            if (riemann_solver_flag .eq. 0) then
               fyro_l = my_ly(i, j, k)
               fymx_l = vy * mx_ly(i, j, k) - by_ly(i, j, k) * bx_ly(i, j, k)
               fymy_l = vy * my_ly(i, j, k) + pr + b2 / 2.0d0 - by_ly(i, j, k) * by_ly(i, j, k)
               fymz_l = vy * mz_ly(i, j, k) - by_ly(i, j, k) * bz_ly(i, j, k)
               fybx_l = vy * bx_ly(i, j, k) - by_ly(i, j, k) * vx
               fyby_l = po_ly(i, j, k)
               fybz_l = vy * bz_ly(i, j, k) - by_ly(i, j, k) * vz
               fyen_l = (en_ly(i, j, k) + pr + b2 / 2.0d0) * vy -                                                                           &
                  (bx_ly(i, j, k) * vx + by_ly(i, j, k) * vy + bz_ly(i, j, k) * vz) * by_ly(i, j, k)
            endif
            fypo_l = ch2 * by_ly(i, j, k)

            vx = mx_ry(i, jp1, k) / ro_ry(i, jp1, k)
            vy = my_ry(i, jp1, k) / ro_ry(i, jp1, k)
            vz = mz_ry(i, jp1, k) / ro_ry(i, jp1, k)
            b2 = bx_ry(i, jp1, k) * bx_ry(i, jp1, k) + by_ry(i, jp1, k) * by_ry(i, jp1, k) + bz_ry(i, jp1, k) * bz_ry(i, jp1, k)
            v2 = vx * vx + vy * vy + vz * vz
            pr = (en_ry(i, jp1, k) - v2 * ro_ry(i, jp1, k) / 2.0d0 - b2 / 2.0d0) * gmm1
            c2 = gm * pr
            s2 = c2 + b2
            ca2 = by_lz(i, jp1, k) * by_lz(i, jp1, k)
            cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_lz(i, jp1, k) / 2.0d0)
            max_speed = max(cfy + abs(vy), max_speed)

            if (riemann_solver_flag .eq. 0) then
               fyro_r = my_ry(i, jp1, k)
               fymx_r = vy * mx_ry(i, jp1, k) - by_ry(i, jp1, k) * bx_ry(i, jp1, k)
               fymy_r = vy * my_ry(i, jp1, k) + pr + b2 / 2.0d0 - by_ry(i, jp1, k) * by_ry(i, jp1, k)
               fymz_r = vy * mz_ry(i, jp1, k) - by_ry(i, jp1, k) * bz_ry(i, jp1, k)
               fybx_r = vy * bx_ry(i, jp1, k) - by_ry(i, jp1, k) * vx
               fyby_r = po_ry(i, j, k)
               fybz_r = vy * bz_ry(i, jp1, k) - by_ry(i, jp1, k) * vz
               fyen_r = (en_ry(i, jp1, k) + pr + b2 / 2.0d0) * vy -                                                                &
                  (bx_ry(i, jp1, k) * vx + by_ry(i, jp1, k) * vy + bz_ry(i, jp1, k) * vz) * by_ry(i, jp1, k)
            endif
            fypo_r = ch2 * by_ry(i, jp1, k)

            if (riemann_solver_flag .eq. 0) then
               fyro(i, j, k) = 0.5d0 * (fyro_l + fyro_r - max_speed * (ro_ry(i, jp1, k) - ro_ly(i, j, k)))
               fymx(i, j, k) = 0.5d0 * (fymx_l + fymx_r - max_speed * (mx_ry(i, jp1, k) - mx_ly(i, j, k)))
               fymy(i, j, k) = 0.5d0 * (fymy_l + fymy_r - max_speed * (my_ry(i, jp1, k) - my_ly(i, j, k)))
               fymz(i, j, k) = 0.5d0 * (fymz_l + fymz_r - max_speed * (mz_ry(i, jp1, k) - mz_ly(i, j, k)))
               fybx(i, j, k) = 0.5d0 * (fybx_l + fybx_r - max_speed * (bx_ry(i, jp1, k) - bx_ly(i, j, k)))
               fyby(i, j, k) = 0.5d0 * (fyby_l + fyby_r - max_speed * (by_ry(i, jp1, k) - by_ly(i, j, k)))
               fybz(i, j, k) = 0.5d0 * (fybz_l + fybz_r - max_speed * (bz_ry(i, jp1, k) - bz_ly(i, j, k)))
               fyen(i, j, k) = 0.5d0 * (fyen_l + fyen_r - max_speed * (en_ry(i, jp1, k) - en_ly(i, j, k)))
            endif
            fypo(i, j, k) = 0.5d0 * (fypo_l + fypo_r - max_speed * (po_ry(i, jp1, k) - po_ly(i, j, k)))

            vx = mx_lz(i, j, k) / ro_lz(i, j, k)
            vy = my_lz(i, j, k) / ro_lz(i, j, k)
            vz = mz_lz(i, j, k) / ro_lz(i, j, k)
            b2 = bx_lz(i, j, k) * bx_lz(i, j, k) + by_lz(i, j, k) * by_lz(i, j, k) + bz_lz(i, j, k) * bz_lz(i, j, k)
            v2 = vx * vx + vy * vy + vz * vz
            pr = (en_lz(i, j, k) - v2 * ro_lz(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
            c2 = gm * pr
            s2 = c2 + b2
            ca2 = bz_lz(i, j, k) * bz_lz(i, j, k)
            cfz = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_lz(i, j, k) / 2.0d0)
            max_speed = cfz + abs(vz)

            if (riemann_solver_flag .eq. 0) then
               fzro_l = mz_lz(i, j, k)
               fzmx_l = vz * mx_lz(i, j, k) - bz_lz(i, j, k) * bx_lz(i, j, k)
               fzmy_l = vz * my_lz(i, j, k) - bz_lz(i, j, k) * by_lz(i, j, k)
               fzmz_l = vz * mz_lz(i, j, k) + pr + b2 / 2.0d0 - bz_lz(i, j, k) * bz_lz(i, j, k)
               fzbx_l = vz * bx_lz(i, j, k) - bz_lz(i, j, k) * vx
               fzby_l = vz * by_lz(i, j, k) - bz_lz(i, j, k) * vy
               fzbz_l = po_lz(i, j, k)
               fzen_l = (en_lz(i, j, k) + pr + b2 / 2.0d0) * vz -                                                                           &
                  (bx_lz(i, j, k) * vx + by_lz(i, j, k) * vy + bz_lz(i, j, k) * vz) * bz_lz(i, j, k)
            endif
            fzpo_l = ch2 * by_lz(i, j, k)

            vx = mx_rz(i, j, kp1) / ro_rz(i, j, kp1)
            vy = my_rz(i, j, kp1) / ro_rz(i, j, kp1)
            vz = mz_rz(i, j, kp1) / ro_rz(i, j, kp1)
            b2 = bx_rz(i, j, kp1) * bx_rz(i, j, kp1) + by_rz(i, j, kp1) * by_rz(i, j, kp1) + bz_rz(i, j, kp1) * bz_rz(i, j, kp1)
            v2 = vx * vx + vy * vy + vz * vz
            pr = (en_rz(i, j, kp1) - v2 * ro_rz(i, j, kp1) / 2.0d0 - b2 / 2.0d0) * gmm1
            c2 = gm * pr
            s2 = c2 + b2
            ca2 = bz_lz(i, j, kp1) * bz_lz(i, j, kp1)
            cfz = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_lz(i, j, kp1) / 2.0d0)
            max_speed = max(cfz + abs(vz), max_speed)

            if (riemann_solver_flag .eq. 0) then
               fzro_r = mz_rz(i, j, kp1)
               fzmx_r = vz * mx_rz(i, j, kp1) - bz_rz(i, j, kp1) * bx_rz(i, j, kp1)
               fzmy_r = vz * my_rz(i, j, kp1) - bz_rz(i, j, kp1) * by_rz(i, j, kp1)
               fzmz_r = vz * mz_rz(i, j, kp1) + pr + b2 / 2.0d0 - bz_rz(i, j, kp1) * bz_rz(i, j, kp1)
               fzbx_r = vz * bx_rz(i, j, kp1) - bz_rz(i, j, kp1) * vx
               fzby_r = vz * by_rz(i, j, kp1) - bz_rz(i, j, kp1) * vy
               fzbz_r = po_rz(i, j, kp1)
               fzen_r = (en_rz(i, j, kp1) + pr + b2 / 2.0d0) * vz -                                                                           &
                  (bx_rz(i, j, kp1) * vx + by_rz(i, j, kp1) * vy + bz_rz(i, j, kp1) * vz) * bz_rz(i, j, kp1)
            endif
            fzpo_r = ch2 * bz_rz(i, j, kp1)

            if (riemann_solver_flag .eq. 0) then
               fzro(i, j, k) = 0.5d0 * (fzro_l + fzro_r - max_speed * (ro_rz(i, j, kp1) - ro_lz(i, j, k)))
               fzmx(i, j, k) = 0.5d0 * (fzmx_l + fzmx_r - max_speed * (mx_rz(i, j, kp1) - mx_lz(i, j, k)))
               fzmy(i, j, k) = 0.5d0 * (fzmy_l + fzmy_r - max_speed * (my_rz(i, j, kp1) - my_lz(i, j, k)))
               fzmz(i, j, k) = 0.5d0 * (fzmz_l + fzmz_r - max_speed * (mz_rz(i, j, kp1) - mz_lz(i, j, k)))
               fzbx(i, j, k) = 0.5d0 * (fzbx_l + fzbx_r - max_speed * (bx_rz(i, j, kp1) - bx_lz(i, j, k)))
               fzby(i, j, k) = 0.5d0 * (fzby_l + fzby_r - max_speed * (by_rz(i, j, kp1) - by_lz(i, j, k)))
               fzbz(i, j, k) = 0.5d0 * (fzbz_l + fzbz_r - max_speed * (bz_rz(i, j, kp1) - bz_lz(i, j, k)))
               fzen(i, j, k) = 0.5d0 * (fzen_l + fzen_r - max_speed * (en_rz(i, j, kp1) - en_lz(i, j, k)))
            endif
            fzpo(i, j, k) = 0.5d0 * (fzpo_l + fzpo_r - max_speed * (po_rz(i, j, kp1) - po_lz(i, j, kp1)))
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  final result
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 2, nz
      km1 = k - 1
      do j = 2, ny
         jm1 = j - 1
         do i = 2, nx
            im1 = i - 1
            ro(i, j, k) = ro(i, j, k) - dt_dx * (fxro(i, j, k) - fxro(im1, j, k)) -                                                &
                                        dt_dy * (fyro(i, j, k) - fyro(i, jm1, k)) -                                                &
                                        dt_dz * (fzro(i, j, k) - fzro(i, j, km1))
            mx(i, j, k) = mx(i, j, k) - dt_dx * (fxmx(i, j, k) - fxmx(im1, j, k)) -                                                &
                                        dt_dy * (fymx(i, j, k) - fymx(i, jm1, k)) -                                                &
                                        dt_dz * (fzmx(i, j, k) - fzmx(i, j, km1))
            my(i, j, k) = my(i, j, k) - dt_dx * (fxmy(i, j, k) - fxmy(im1, j, k)) -                                                &
                                        dt_dy * (fymy(i, j, k) - fymy(i, jm1, k)) -                                                &
                                        dt_dz * (fzmy(i, j, k) - fzmy(i, j, km1))
            mz(i, j, k) = mz(i, j, k) - dt_dx * (fxmz(i, j, k) - fxmz(im1, j, k)) -                                                &
                                        dt_dy * (fymz(i, j, k) - fymz(i, jm1, k)) -                                                &
                                        dt_dz * (fzmz(i, j, k) - fzmz(i, j, km1))
            bx(i, j, k) = bx(i, j, k) - dt_dx * (fxbx(i, j, k) - fxbx(im1, j, k)) -                                                &
                                        dt_dy * (fybx(i, j, k) - fybx(i, jm1, k)) -                                                &
                                        dt_dz * (fzbx(i, j, k) - fzbx(i, j, km1))
            by(i, j, k) = by(i, j, k) - dt_dx * (fxby(i, j, k) - fxby(im1, j, k)) -                                                &
                                        dt_dy * (fyby(i, j, k) - fyby(i, jm1, k)) -                                                &
                                        dt_dz * (fzby(i, j, k) - fzby(i, j, km1))
            bz(i, j, k) = bz(i, j, k) - dt_dx * (fxbz(i, j, k) - fxbz(im1, j, k)) -                                                &
                                        dt_dy * (fybz(i, j, k) - fybz(i, jm1, k)) -                                                &
                                        dt_dz * (fzbz(i, j, k) - fzbz(i, j, km1))
            en(i, j, k) = en(i, j, k) - dt_dx * (fxen(i, j, k) - fxen(im1, j, k)) -                                                &
                                        dt_dy * (fyen(i, j, k) - fyen(i, jm1, k)) -                                                &
                                        dt_dz * (fzen(i, j, k) - fzen(i, j, km1))
            po(i, j, k) = po(i, j, k) - dt_dx * (fxpo(i, j, k) - fxpo(im1, j, k)) -                                                &
                                        dt_dy * (fypo(i, j, k) - fypo(i, jm1, k)) -                                                &
                                        dt_dz * (fzpo(i, j, k) - fzpo(i, j, km1))
         enddo
      enddo
   enddo

   call resistivity(bx_resis, by_resis, bz_resis, en_resis,                                                                        &
      ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, x, y, z, dx, dy, dz, nx, ny, nz)

   call source_term (mx, my, mz, bx, by, bz, en, po,                                                                               &
      ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, po_pre, gx, gy, gz, bx_resis, by_resis, bz_resis, en_resis,          &
      ch2, cp2, dx, dy, dz, dt, nx, ny, nz)

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
      fxro_bnd(j, k, 1) = fxro(ng, j, k)
      fxro_bnd(j, k, 2) = fxro(nx - ng, j, k)

      fxmx_bnd(j, k, 1) = fxmx(ng, j, k)
      fxmx_bnd(j, k, 2) = fxmx(nx - ng, j, k)

      fxmy_bnd(j, k, 1) = fxmy(ng, j, k)
      fxmy_bnd(j, k, 2) = fxmy(nx - ng, j, k)

      fxmz_bnd(j, k, 1) = fxmz(ng, j, k)
      fxmz_bnd(j, k, 2) = fxmz(nx - ng, j, k)

      fxbx_bnd(j, k, 1) = fxbx(ng, j, k)
      fxbx_bnd(j, k, 2) = fxbx(nx - ng, j, k)

      fxby_bnd(j, k, 1) = fxby(ng, j, k)
      fxby_bnd(j, k, 2) = fxby(nx - ng, j, k)

      fxbz_bnd(j, k, 1) = fxbz(ng, j, k)
      fxbz_bnd(j, k, 2) = fxbz(nx - ng, j, k)

      fxen_bnd(j, k, 1) = fxen(ng, j, k)
      fxen_bnd(j, k, 2) = fxen(nx - ng, j, k)

      fxpo_bnd(j, k, 1) = fxpo(ng, j, k)
      fxpo_bnd(j, k, 2) = fxpo(nx - ng, j, k)
   enddo
   enddo

   do k = 1, nz
   do i = 1, nx
      fyro_bnd(i, k, 1) = fyro(i, ng, k)
      fyro_bnd(i, k, 2) = fyro(i, ny - ng, k)

      fymx_bnd(i, k, 1) = fymx(i, ng, k)
      fymx_bnd(i, k, 2) = fymx(i, ny - ng, k)

      fymy_bnd(i, k, 1) = fymy(i, ng, k)
      fymy_bnd(i, k, 2) = fymy(i, ny - ng, k)

      fymz_bnd(i, k, 1) = fymz(i, ng, k)
      fymz_bnd(i, k, 2) = fymz(i, ny - ng, k)

      fybx_bnd(i, k, 1) = fybx(i, ng, k)
      fybx_bnd(i, k, 2) = fybx(i, ny - ng, k)

      fyby_bnd(i, k, 1) = fyby(i, ng, k)
      fyby_bnd(i, k, 2) = fyby(i, ny - ng, k)

      fybz_bnd(i, k, 1) = fybz(i, ng, k)
      fybz_bnd(i, k, 2) = fybz(i, ny - ng, k)

      fyen_bnd(i, k, 1) = fyen(i, ng, k)
      fyen_bnd(i, k, 2) = fyen(i, ny - ng, k)

      fypo_bnd(i, k, 1) = fypo(i, ng, k)
      fypo_bnd(i, k, 2) = fypo(i, ny - ng, k)
   enddo
   enddo

   do j = 1, ny
   do i = 1, nx
      fzro_bnd(i, j, 1) = fzro(i, j, ng)
      fzro_bnd(i, j, 2) = fzro(i, j, nz - ng)

      fzmx_bnd(i, j, 1) = fzmx(i, j, ng)
      fzmx_bnd(i, j, 2) = fzmx(i, j, nz - ng)

      fzmy_bnd(i, j, 1) = fzmy(i, j, ng)
      fzmy_bnd(i, j, 2) = fzmy(i, j, nz - ng)

      fzmz_bnd(i, j, 1) = fzmz(i, j, ng)
      fzmz_bnd(i, j, 2) = fzmz(i, j, nz - ng)

      fzbx_bnd(i, j, 1) = fzbx(i, j, ng)
      fzbx_bnd(i, j, 2) = fzbx(i, j, nz - ng)

      fzby_bnd(i, j, 1) = fzby(i, j, ng)
      fzby_bnd(i, j, 2) = fzby(i, j, nz - ng)

      fzbz_bnd(i, j, 1) = fzbz(i, j, ng)
      fzbz_bnd(i, j, 2) = fzbz(i, j, nz - ng)

      fzen_bnd(i, j, 1) = fzen(i, j, ng)
      fzen_bnd(i, j, 2) = fzen(i, j, nz - ng)

      fzpo_bnd(i, j, 1) = fzpo(i, j, ng)
      fzpo_bnd(i, j, 2) = fzpo(i, j, nz - ng)
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine tvd_lf
